An Average Flow Model of the Reynolds Roughness Including a Mass-Flow Preserving Cavitation Model

An average Reynolds equation for predicting the effects of deterministic periodic roughness, taking JFO mass flow preserving cavitation model into account, is introduced based upon double scale analysis approach. This average Reynolds equation can be used both for a microscopic interasperity cavitation and a macroscopic one. The validity of such a model is verified by numerical experiments both for one dimensional and two dimensional roughness patterns. 1 NOMENCLATURE A", B", Ai, Bi = partial differential operators a?ij , a?i , a0i = auxiliary homogenized coefficients A?ij , B? i , B0 i = homogenized coefficients h1, h2 = description of the gap h, h" = actual gap hs = smooth part of the gap hr = amplitude of the roughness p = pressure p0, p1... = approximations of the pressure Q = input flow value U = velocity x = (x1; x2) = dimensionless space coordinates y = (y1; y2) = microscale coordinates X = (X1; X2) = oblique coordinates X 0 = (X 0 1; X 0 2) = real coordinates Y =℄0; 1[ ℄0; 1[ = rescaled microcell = obliqueness angle = n = normal derivative = viscosity " = roughness spacing = saturation 0 = microscopic homogenized saturation , 1, 2 = macrohomogenized saturations wi, 0i = auxiliary functions defined on Y Y = average operator with respect to y [ ℄Y1 = average operator with respect to y1 [ ℄Y2 = average operator with respect to y2 0 Introduction The effects of the surface roughness on the behavior of a thin film flow has long been the subject of intensive studies. Various ways have been introduced to study Reynolds roughness by seeking an average equation with smooth coefficients. Some of the most popular results are the Christensen formula [1] for longitudinal and transverse roughness and the Patir and Cheng flow factor model [2] for a more general surface roughness pattern. Two wide classes of results can be outlined. In the first one, which is deterministic, a periodic description of the surfaces is often assumed to be known and linked to a specific process of the surface [3]. It is possible to distinguish macrovariables and microvariables and to use a mathematical homogenization approach to rigorously obtain an average Reynolds equation by making the period of the roughness tend to zero [4]. The coefficients of this average Reynolds equation implicitly contain the description of the microroughness elementary cell. The second class of results deals with a statistical description of the surface roughness. Following the Patir and Cheng approach, numerous authors proposed an average Reynolds equation in which the coefficients included the knowledge of the surface statistics by way of flow factors which can be evaluated by numerical experiments. Rigorously speaking, this approach is less satisfactory than the first one, assuming a priori the existence of a control volume in which 2 the average flow rates can be equivalently expressed in terms of flow factors. The number and quantities (Peklenik number, combined root mean square roughness...) involved in the characterization of the flow factors can also be discussed. Moreover, as the initial Reynolds equation, the average Reynolds equation can be expressed in terms of r (Krp) = F in which K is a diagonal matrix. This seems to be contradictory with the result obtained by the first approach in which K is a non diagonal matrix for two dimensional general roughness pattern [5]. Up to now, these averaging processes never take cavitation into account. A common procedure is to use the average equation instead of the classical Reynolds equation with Gumbel and Swift Steiber boundary conditions or to include it in the S.O.R. algorithm proposed by Richardson, thus obtaining the splitting of the lubricated device in two areas. In a first area, the pressure is greater than the cavitation pressure and the average Reynolds equation is valid; in the other area, pressure is equal to the cavitation pressure. It is well known [6–8], however, that none of these models is mass preserving, especially through the cavitation area. Jakobsson, Floberg and Olsson (JFO) [9,10] developed a set of conditions for the cavitation boundary that properly takes the conservation of mass into account in the entire device. Elrod [11, 12] proposed a slightly modified formulation and a related specific algorithm. The mathematical related problem evidence a hyperbolic-parabolic feature which renders difficult both theoretical study and numerical experiments [7, 13–15]. It is the goal of this paper to develop in a rigorous way an average JFO Reynolds equation for the deterministic periodic roughness pattern. So far, few papers have been devoted to such a problem. Recently, the interasperity cavitation has been studied by way of a statistical approach [16, 17]. The Patir and Cheng flow factor method is extended and an average Reynolds equation is proposed. The resulting equation has the same left-hand side that in the Patir and Cheng equation (cavitation has no effect on the corresponding flow factors) while the right-hand side of the equation is modified and new flow factors are introduced. At last Harp and Salant [17] proposed to modify the boundary conditions by a value which is a function of the wavelength of the roughness. Our approach is quite different and explicitly based upon the introduction of fast and slow variables. The initial equation is rewritten in terms of these two variables and asymptotic expansion of the pressure is introduced with respect to a small parameter associated to the roughness wavelength. The goal is to find an equation satisfied by the first terms of the expansion. Some assumptions about the shape of the roughness appear to be necessary to solve the problem, leading to a new average Reynolds cavitation equation. This equation has numerous common features with the initial Reynolds equation: it is also a two unknowns pressure-saturation formulation. Some particular cases transverse, longitudinal roughness patterns will be studied in details. 1 Basic equations Our studied cavitation model, like the Elrod algorithm and its variants, views the film as a mixture. It does not, however, make the assumption of liquid compressibility in the full film area as in [15] and some other papers. As in [18, 19], only the liquid-vapor mixture in the cavitated region is assumed compressible. The flow obeys the following “universal” Reynolds equation (here written in a dimensionless form) through all the gap in which the pressure cavitation is assumed to be zero in the cavitation area 2 Xi=1 xi h3 p xi! = h x1 ; (1) p 0; (2) 0 1; (3) 3 p (1 ) = 0: (4) In this steady state isoviscous version of the equation, p is the pressure, is the relative mixture density, h the film thickness, x1 is the direction of the effective relative velocity of the shaft, while x2 is the transverse direction. This system of equations can be understood as follows (see [7,9–12,14,18] for various comments and meaning of the variable): the well-known Reynolds equation holds in the full film region, that is p > 0 and = 1, a mass flow conserving equation h= x1 = 0 holds in the cavitated region with p = 0 and 0 < < 1. a boundary condition which is also mass flow preserving at the (unknown) interface between the two regions: h3 p n + h os(n; x1) = h os(n; x1): The reason to retain this specific cavitation equation is that it has been the subject of numerous mathematical studies [7] giving a strong and rigorous basis to the following manipulations [20]. To be noticed, however, that our approach can be applied without difficulty to other cavitation models as the one in [15]. Last, it has to be mentioned that this equation takes both macrocavitation (associated to the occurrence of a diverging part of a bearing for example) and interasperity cavitation into account. The boundary conditions depend on the considered device. However, the following ones are often used, corresponding for example to a journal bearing with an axial supply groove. The pressure is imposed at two circumferential locations and one axial location. The last boundary condition is an input flow condition at the axial location corresponding to the supply groove: (x)h(x) h3(x) p x1 (x) = Q: (5) For small values of Q, starvation may occur in the vicinity of the supply groove. 2 Asymptotic expansion Let us suppose that the roughness is periodically reproduced in the two x1 and x2 directions from an elementary cell Y (or “miniature bearing” in Tonder’s terminology). We denote by " the ratio of the homothetic transformation passing from the elementary cell Y = Y1 Y2 to the real bearing and by y1 = x1=" and y2 = x2=" the local variables (see FIG. 1). Let us now consider shapes that can be written as h"(x) = h(x; x="). We suppose furthermore that they are described as h"(x) = h1 x; x1 " h2 x; x2 " which allows us to take into account either transverse or longitudinal roughness, but also more general two dimensional roughness. Introducing now the fast variables y1 and y2, it appears that the new expression for the gap 4 "Y Y y1 y2 x1 x2 y = x" Figure 1. Macroscopic domain and elementary cells is: h(x; y) = h1 (x; y1)h2 (x; y2) : (6) The combined computation in terms of (x1; x2) or (y1; y2) is an important feature of the method. It is convenient to consider first x and y as independent variables and to replace next y by x=" (see [4]). 2.1 Formulation of average equations We denote by A" the initial differential Reynolds operator A"[ ℄ = 2 X j=1 xj h3 x; x" [ ℄ xj ! ; and we also define the right-hand side operator B"[ ℄ = x1 h x; x" [ ℄ : The Reynolds equation (1) becomes A"(p) = B"( ): 5 The underscript " indicates the dependance of the real pressure on the microtexture related to ". We also define the following operators: A1[ ℄ = 2 X j=1 yj h3 (x; y) [ ℄ yj ! ; A2[ ℄ = 2 X j=1 yj h3(x; y) [ ℄ xj !+ 2 X j=1 xj h3(x; y) [ ℄ yj ! ; A3[ ℄ = 2 X j=1 xj h3(x; y) [ ℄ xj ! ; and also Bi 1[ ℄ = yi (h(x; y) [ ℄ ) ; i = 1; 2; Bi 2[ ℄ = xi (h(x; y) [ ℄ ) ; i = 1; 2: If applied to a function of (x; x="), the operators become A" = 1="2 A1 + 1=" A2 + A3 ; (7) B" = 1=" B1 1 +B1 2 : (8) We shall look for an asymptotic expansion of the solutions p(x) = p0(x; x" ) + "p1 x; x" + "2p2 x; x" + :::; (9) (x) = 0 x; x" ; (10) each unknown pi and 0 being a function of (x; y). The problem of the boundary conditions to be satisfied by the pi is somewhat difficult but may be summarized as follows. (i) The natural boundary conditions on (p"; ") are assigned to p0 and an equivalent saturation linked to 0, which will be developped in next subsection. (ii) The function pi, i 1, are Y periodic, i.e. periodic in the two variables y1, y2, for each value of (x1; x2). To be noticed that unlike of p, we do not introduce an asymptotic expansion for . This can be explained by observing the evolution of p and as " tends to 0 (see FIG.2 for instance). Clearly, the oscillations of the pressure are decreasing and p tends to a smooth function (namely p0 which, actually, does not depend on the fast variable as it will be pointed out further). This is not the case for and an asymptotic smooth limit cannot be considered. 6 We shall see later that the functions pi, i 1, are defined up to an additive constant. Moreover, from Equations (2)–4), the following properties hold: p0(x; y) 0; (11) 0 0(x; y) 1; (12) p0(x; y) (1 0(x; y)) = 0: (13) Putting Equations (9) and (10) into Equation (1) and taking account of Equations (7) and (8), one can write by an identification procedure: A1p0 = 0; (14) A1p1 + A2p0 = B1 1 0; (15) A1p2 + A2p1 + A3p0 = B1 2 0: (16) Let us remark that these equations are of the following type: For a given F , find a function q, depending on the variable y, q being Y periodic, such that (x is a parameter), A1q = F: (17) A condition to have a solution for Equation (17) is ZY F (x; y)dy = 0: (18) Moreover, if q is a solution, then q+ with any constant with respect to y is also a solution. Applying Condition (18) to Equation (14), we deduce that p0 does not depend on y p0(x): (19) Let us suppose now that p0 is known, and noticing that, due to boundary conditions, (B1 1 0 A2p0) satisfies Equation (18), existence of p1 is guaranteed. Now we can represent p1 as a function of p0 in a more usable form. We define wi and 0i (i = 1; 2) as the Y periodic solutions (up to an additive constant) of the following local problems: A1 wi = h3 yi ; i = 1; 2; (20) A1 0i = 0h yi ; i = 1; 2: (21) 7 The solution of Equation (15) reduces to p1(x; y) = 01(x; y) p0 x1 (x)w1(x; y) p0 x2 (x)w2(x; y): (22) The same procedure can be used to ensure the existence of p2, but in that step, the corresponding condition (18) applied to Equation (16) becomes ZY (B1 2 0 A2p1 A3p0) dy = 0: (23) Then the main idea is to put Equation (22) into Equation (23), so that the only remaining unknowns are p0 and 0. By analogy with the probabilistic framework, we denote by uY the local average of any Y periodic function u: uY (x) = 1 [Y ℄ ZY u(x; y) dy: By exchanging the integral and the derivation symbols, and after some calculations, Equation (23) becomes Xi;j xi A?ij p0 xj! = B0 1 x1 + B0 2 x2 ! ; (24) where (i; j = 1; 2 and j 6 = i) A?ii = h3Y h3 wi yi Y ; A?ij = h3 wj yi Y = h3 wi yj Y = A?ji; and also B0 1 = 0hY h3 01 y1 Y ; B0 2 = h3 01 y2 Y : Equation (24) deals with any periodic roughness pattern. To be noticed is the fact that the differential operator is no more of the Reynolds type since extra terms 2p0= xi xj appear. The right-hand side also contains an additive term in the x2 direction. However, the link between p0 and 0 is not so clear. This is a major obstacle which prevents from getting a tractable equation. Nevertheless, Assumption 6 allows us to solve the following difficulties: 8 Computation of A?ii, i = 1; 2: Let us recall Equation (20) with i = 1: y1 h3 w1 y1 !+ y2 h3 w1 y2 ! = h3 y1 : Since h3 w1= y2 is Y periodic, averaging this equation over Y2 gives y1 "h3 w1 y1 #Y2! = [h3℄Y2 y1 ; where [ ℄[Yi℄ is the averaging operator over Yi (for i = 1; 2). Thus we have, by integrating in the y1 variable and using Equations(6): "h3 h3 w1 y1 #Y2 = C1; where C1 is a constant with respect to y. Let us notice that, averaging the earlier equation over Y1 simply gives C1 = A?11. Thus, it remains to calculate C1. Dividing each side of the previous equation by h31: hh32iY2 "h32 w1 y1 #Y2 = C1 h31 and, since w1 is Y periodic, averaging over Y1 gives A?11 = h32Y h 3 1 Y : (25) Following the same procedure, we state: A?22 = h31Y h 3 2 Y : (26) Computation of A?ij , i 6= j: Starting from Equation (20) with i = 1, since h3 h3 w1= y1 is Y periodic, averaging this equation over Y1 gives y2 "h3 w1 y2 #Y1! = 0: 9 Thus we have, by integrating in the y2 variable "h3 w1 y2 #Y1 = C2; where C2 is a constant with respect of y. Similarly to the computation of A?ii, one has C2 = A?12 = A?21. Dividing each side of the equation by h32: C2 h32 = "h31 w1 y2 #Y1 ; and, since w1 is Y periodic, averaging over Y2 gives C2a 1 1 Y = 0, i.e. A?12 = A?21 = 0: (27) Now, it remains to calculate the right-hand side of the Reynolds equation. Computation of B0 1 : Let us recall Equation (21) with i = 1: y1 h3 01 y1 !+ y2 h3 01 y2 ! = 0h y1 : Since h3 01= y2 is Y periodic, averaging this equation over Y2 gives y1 "h3 01 y1 #Y2! = [ 0h℄Y2 y1 : Thus we have, by integrating in the y1 variable: " 0h h3 10 y1 #Y2 = C3; where C3 is a constant with respect to y. Clearly, we have C3 = B0 1 . Dividing each side of the equation by h31: " 0h h31 #Y2 "h32 01 y1 #Y2 = C3 h31 ; 10 and, since 01 is Y periodic, averaging over Y1 gives 0h=h31Y = C3h 3 1 Y , i.e. B0 1 = 0h2 h21 !Y h 3 1 Y (28) Computation of B0 2 : Starting from Equation (21) with i = 1, since the function h3 h3 01= y1 is Y periodic, averaging this equation over Y1 gives y2 "h3 01 y2 #Y1! = 0: Thus we have, by integrating in the y2 variable: "h3 01 y2 #Y1 = C4; where C4 is a constant with respect of y. We have C = B0 2 . Then dividing each side by h32: C4 h32 = "h31 01 y2 #Y1 ; and, since 01 is Y periodic, averaging over Y2 gives C4h 3 1 Y = 0, i.e. B0 2 = 0: (29) Now, it is obvious that Equation (24) can be written in a more simple way by using Equations (25)–(29). Before that, let us write the term B0 1 in a more usable form. Defining the quantities B? 1 = h 2 1 Y h 3 1 Y h2Y ; (30) = 1 h2Y h 2 1 Y 0h2 h21 !Y ; (31) we get B0 1 = B? 1 . Moreover, from Equations (12) and (13), we immediately have: 0 (x) 1; (32) p0(x) (1 (x)) = 0; (33) 11 so that the homogenized equations appear to be 2 Xi=1 xi A?ii p0 xi! = B? 1 x1 ; (34) p0 0; (35) 0 1; (36) p0 (1 ) = 0; (37) where A?11, A?22 and B? 1 are, respectively, given by Equations (25), (26) and (30). Moreover, the link between a new (smooth) “macroscopic” saturation and the (oscillating) “microscopic” saturation 0 is given by Equation (31). As an important feature, is not the average of the microscopic saturation 0. 2.2 Average boundary condition When the pressure is imposed, the corresponding average boundary condition is assigned to p0. When an input flow is given on a supply line, the average flow condition is obtained following the asymptotic expansion method. Taking account of roughness patterns, Equation (5) becomes: (x)h x; x" h3 x; x" p x1 (x) = Q: (38) Putting Equations (9) and (10) into Equation (38), one can write by an identification procedure: 0(x; y)h(x; y) h3(x; y) p0 x1 (x) + p1 y1 (x; y)! = Q: Putting Equation (22) into it gives 0h h3 01 y1 ! h3 h3 w1 y1 ! p0 x1 + h3 w2 y1 ! p0 x2 = Q: Averaging over Y gives the boundary condition relating p0 and at the supply groove: B0 1 A?11 p0 x1 A?12 p0 x2 = Q; and since A?12 = 0 and B0 1 = B? 1 , one gets: B? 1 A?11 p0 x1 = Q: (39) The next subsection deals with two main particular cases: transverse or longitudinal roughness. 12 2.3 Particular cases Transverse roughness: when the roughness does not depend on y2, we have the homogenized equation, easily deduced from Equations (25)–(31) x1 1 h 3Y p0 x1!+ x2 h3Y p0 x2! = x1 0 h 2Y h 3Y 1A ; with = 1 h 2Y 0 h2!Y and the boundary condition at the supply groove, deduced from Equation (39), should be read as: h 2Y h 3Y 1 h 3Y p0 x1 = Q: Longitudinal roughness: when the roughness does not depend on y1, we get x1 h3Y p0 x1!+ x2 1 h 3Y p0 x2! = x1 hY ; with = 0hY hY , and the boundary condition at the supply groove should be read as: hY h3Y p0 x1 = Q: 3 Oblique roughness Let us consider gaps that can be written as: h"(x) = h1 x; X1(x) " !h2 x; X2(x) " ! ; with (X1(x) = os x1 + sin x2; X2(x) = sin x1 + os x2; which allows us to take into account oblique roughness (with h2 1 for instance). The idea is to introduce a change of coordinates so that the assumption of Section B.2 on the roughness form in the new coordinates system is valid. The first step is to rewrite Equation (1) in the X coordinates: 2 Xi=1 Xi h3" p Xi! = h" X1 os h" X2 sin ! : 13 Working now in the X coordinates and using the operators defined in Section B.2 (up to the writing in the X coordinates), we apply the asymptotic expansion technique to the earlier equation. With the formal asymptotic expansion used in Section B.2, we have in the (X; y) coordinates (with y = X="): A1p0 = 0; (40) A1p1 + A2p0 = B1 1 0 os B2 1 0 sin ; (41) A1p2 + A2p1 + A3p0 = B1 2 0 os B2 2 0 sin : (42) As in Section B.2, p0 only depends on the X variable. Equation (41) allows us to determine p1: p1(X; y) = 01(X; y) os 02(X; y) sin w1(X; y) p0 X1 (X) w2(X; y) p0 X2 (X): Then, putting the earlier expression into Equation (42) gives: Xi;j Xi a?ij p0 Xj! = X1 b011 os + b012 sin + X2 b021 os + b022 sin ; where the coefficients, which are easily computed as in Section B.2, are given by (i; j = 1; 2, i 6= j): a?ii = h3Y h3 wi yi Y = h3j Y h 3 i Y ; (43) a?ij = h3 wj yi Y = 0; (44) and also (i; j = 1; 2, i 6= j) b0ii = 0hY h3 0i yi Y = 1 h 3 i Y 0hj h2i !Y ; (45) b0ij = h3 0i yj Y = 0: (46) Finally, as in Section B.2, defining the quantities b?i = h 2 i Y h 3 i Y hjY ; (47) i = 1 hjY h 2 i Y 0hj h2i !Y ; (48)